GH_2gen_2021
Mohammad Inam Jameel
4/12/2022
Relative contributions
SLA individual for full dataset
RC_SLA<- lmer(SLA100 ~ treatment*elev*mat_treat+S_init_size +(1|block),control = lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e7)), data = gh2)
Anova(RC_SLA)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: SLA100
## Chisq Df Pr(>Chisq)
## treatment 13.1003 1 0.0002952 ***
## elev 13.6912 1 0.0002155 ***
## mat_treat 7.3609 2 0.0252121 *
## S_init_size 38.9269 1 4.4e-10 ***
## treatment:elev 2.3215 1 0.1275961
## treatment:mat_treat 1.6117 2 0.4467191
## elev:mat_treat 1.3702 2 0.5040389
## treatment:elev:mat_treat 0.0353 2 0.9825290
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot (RC_SLA)
visreg(RC_SLA, overlay = TRUE, "elev", by="treatment", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="Specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) 
visreg(RC_SLA,overlay = TRUE, "elev", by="mat_treat",type="conditional", scale = "response",
xlab="Elevation (m)", ylab="Specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
#fill=list(col=grey(c(0.7,0.9))),
#line=list(col=grey(c(0,0.5))),
#points=list(cex=1.5,col=grey(c(0.2,0.8)))
) 
SLA family level
#read in emmeans
LSmeans_SLA <- read.csv("LSmeans_SLA_2.csv", stringsAsFactors=TRUE)
RC_SLA_means<- lm(emmean ~ elevation*treatment*mat_treat, data = LSmeans_SLA)
Anova(RC_SLA_means)## Anova Table (Type II tests)
##
## Response: emmean
## Sum Sq Df F value Pr(>F)
## elevation 13411 1 6.7190 0.01095 *
## treatment 73602 1 36.8764 2.229e-08 ***
## mat_treat 5372 2 1.3458 0.26497
## elevation:treatment 2232 1 1.1185 0.29277
## elevation:mat_treat 2575 2 0.6450 0.52680
## treatment:mat_treat 690 2 0.1729 0.84151
## elevation:treatment:mat_treat 22 2 0.0056 0.99441
## Residuals 201587 101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(RC_SLA_means)
##use this
visreg(RC_SLA_means, overlay = TRUE, "elevation", by="treatment", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) 
visreg(RC_SLA_means, overlay = FALSE, "elevation", by="mat_treat", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
#fill=list(col=grey(c(0.7,0.9))),
#line=list(col=grey(c(0,0.5))),
#points=list(cex=1.5,col=grey(c(0.2,0.8)))
) 
Trichomes individual
RC_Tri<- glmer(trichomes ~ mat_treat*elev*treatment+S_init_size +(1|block),family = Gamma(link=log), data = gh2) ## boundary (singular) fit: see ?isSingular#had errors but this seems to work when i specify link=log
#singular fit
Anova(RC_Tri) # elev is sig, sig interaction btw mat_treat and elev## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: trichomes
## Chisq Df Pr(>Chisq)
## mat_treat 9.6468 2 0.008039 **
## elev 111.6048 1 < 2.2e-16 ***
## treatment 2.2309 1 0.135277
## S_init_size 15.4506 1 8.469e-05 ***
## mat_treat:elev 10.5228 2 0.005188 **
## mat_treat:treatment 1.5863 2 0.452418
## elev:treatment 0.1045 1 0.746480
## mat_treat:elev:treatment 0.7294 2 0.694389
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot (RC_Tri)
visreg(RC_Tri, "elev", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="Trichome density", partial=FALSE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) ## Warning: Note that you are attempting to plot a 'main effect' in a model that contains an
## interaction. This is potentially misleading; you may wish to consider using the 'by'
## argument.## Conditions used in construction of plot
## mat_treat: Parent
## treatment: Control
## S_init_size: 0.251658
## block: 9
visreg(RC_Tri, "elev", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="Trichome density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) ## Warning: Note that you are attempting to plot a 'main effect' in a model that contains an
## interaction. This is potentially misleading; you may wish to consider using the 'by'
## argument.## Conditions used in construction of plot
## mat_treat: Parent
## treatment: Control
## S_init_size: 0.251658
## block: 9
visreg(RC_Tri, "elev", by="mat_treat", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="Trichome density", partial=FALSE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) 
visreg(RC_Tri, "elev", by="mat_treat", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="ç", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) 
Trichomes family level
LSmeans_tri <- read.csv("LSmeans_tri.csv", stringsAsFactors=TRUE)
RC_tri_means<- glm(emmean ~ mat_treat*elev*treatment, data = LSmeans_tri)
Anova(RC_tri_means)## Analysis of Deviance Table (Type II tests)
##
## Response: emmean
## LR Chisq Df Pr(>Chisq)
## mat_treat 0.5249 2 0.7691807
## elev 15.0428 1 0.0001051 ***
## treatment 0.0234 1 0.8782974
## mat_treat:elev 2.6707 2 0.2630689
## mat_treat:treatment 0.3957 2 0.8205087
## elev:treatment 0.0217 1 0.8827648
## mat_treat:elev:treatment 0.1180 2 0.9427011
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(RC_tri_means)
##use this
visreg(RC_tri_means, overlay = TRUE, "elev", by="treatment", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="Trichome density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) 
visreg(RC_tri_means, overlay = FALSE, "elev", by="mat_treat", type="conditional", scale = "response",
xlab="Elevation (m)", ylab="Trichome density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
#fill=list(col=grey(c(0.7,0.9))),
#line=list(col=grey(c(0,0.5))),
#points=list(cex=1.5,col=grey(c(0.2,0.8)))
) 
Plasticity overview
SLA
#read in emmeans. using lsmeans because its easier to see the graph
LSmeans_SLA <- read.csv("LSmeans_SLA_2.csv", stringsAsFactors=TRUE)
SLA <- lm(emmean ~ treatment*mat_treat, data = LSmeans_SLA)
Anova (SLA)## Anova Table (Type II tests)
##
## Response: emmean
## Sum Sq Df F value Pr(>F)
## treatment 74246 1 36.1403 2.587e-08 ***
## mat_treat 5233 2 1.2736 0.2840
## treatment:mat_treat 652 2 0.1586 0.8535
## Residuals 219818 107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##use this
visreg(SLA, overlay = TRUE, "mat_treat", by="treatment", type="conditional", scale = "response",
xlab="Maternal environment", ylab="Specific leaf are", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) 
#plasticity figure
ggplot(
#subset(
LSmeans_SLA,
# treatment == "Herbivorized"),
aes(x=treatment, y=emmean, color = genotype)) + # Change fill to color
theme_classic(base_size = 22) +
#geom_point() +
stat_summary(fun=mean, position = "dodge") +
#stat_summary(
# geom="errorbar",
# fun.data= mean_cl_boot,
# width = 0.1, size = 0.2, col = "grey57"
# ) +
# Lines by species using grouping
#scale_color_viridis(discrete = FALSE) +
#scale_fill_manual(values = cbp2)+
stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
ylab("Specific Leaf Area") +
xlab("Treatment")## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Removed 1 rows containing non-finite values (stat_summary).## Warning: Width not defined. Set with `position_dodge(width = ?)`## Warning: Removed 38 rows containing missing values (geom_segment).
#plasticity figure
ggplot(
subset(
LSmeans_SLA,
treatment == "Herbivorized"),
aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
theme_classic(base_size = 22) +
#geom_point() +
stat_summary(fun=mean, position = "dodge") +
#stat_summary(
# geom="errorbar",
# fun.data= mean_cl_boot,
# width = 0.1, size = 0.2, col = "grey57"
# ) +
# Lines by species using grouping
#scale_color_viridis(discrete = FALSE) +
#scale_fill_manual(values = cbp2)+
stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
ylab("Specific Leaf Area") +
xlab("Maternal enviornment")## Warning: Width not defined. Set with `position_dodge(width = ?)`## Warning: Removed 57 rows containing missing values (geom_segment).
#plasticity figure
ggplot(
subset(
LSmeans_SLA,
treatment == "Control"),
aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
theme_classic(base_size = 22) +
#geom_point() +
stat_summary(fun=mean, position = "dodge") +
#stat_summary(
# geom="errorbar",
# fun.data= mean_cl_boot,
# width = 0.1, size = 0.2, col = "grey57"
# ) +
# Lines by species using grouping
#scale_color_viridis(discrete = FALSE) +
#scale_fill_manual(values = cbp2)+
stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
ylab("Specific Leaf Area") +
xlab("Maternal enviornment")## Warning: Removed 1 rows containing non-finite values (stat_summary).## Warning: Removed 1 rows containing non-finite values (stat_summary).## Warning: Width not defined. Set with `position_dodge(width = ?)`## Warning: Removed 56 rows containing missing values (geom_segment).
Trichome Density
#read in emmeans. using lsmeans because its easier to see the graph
LSmeans_tri <- read.csv("LSmeans_tri.csv", stringsAsFactors=TRUE)
TRI <- glm(emmean ~ treatment*mat_treat, data = LSmeans_tri)
Anova (TRI)## Analysis of Deviance Table (Type II tests)
##
## Response: emmean
## LR Chisq Df Pr(>Chisq)
## treatment 0.00821 1 0.9278
## mat_treat 0.43504 2 0.8045
## treatment:mat_treat 0.28404 2 0.8676##use this
visreg(TRI, overlay = TRUE, "mat_treat", by="treatment", type="conditional", scale = "response",
xlab="Maternal environment", ylab="Trichome Density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col=grey(c(0.7,0.9))),
line=list(col=grey(c(0,0.5))),
points=list(cex=1.5,col=grey(c(0.2,0.8)))) 
#plasticity figure
ggplot(
#subset(
LSmeans_tri,
# treatment == "Herbivorized"),
aes(x=treatment, y=emmean, color = genotype)) + # Change fill to color
theme_classic(base_size = 22) +
#geom_point() +
stat_summary(fun=mean, position = "dodge") +
#stat_summary(
# geom="errorbar",
# fun.data= mean_cl_boot,
# width = 0.1, size = 0.2, col = "grey57"
# ) +
# Lines by species using grouping
#scale_color_viridis(discrete = FALSE) +
#scale_fill_manual(values = cbp2)+
stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
ylab("Trichome Density") +
xlab("Treatment")## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Removed 1 rows containing non-finite values (stat_summary).## Warning: Width not defined. Set with `position_dodge(width = ?)`## Warning: Removed 38 rows containing missing values (geom_segment).
#plasticity figure
ggplot(
subset(
LSmeans_tri,
treatment == "Herbivorized"),
aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
theme_classic(base_size = 22) +
#geom_point() +
stat_summary(fun=mean, position = "dodge") +
#stat_summary(
# geom="errorbar",
# fun.data= mean_cl_boot,
# width = 0.1, size = 0.2, col = "grey57"
# ) +
# Lines by species using grouping
#scale_color_viridis(discrete = FALSE) +
#scale_fill_manual(values = cbp2)+
stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
ylab("Trichome Density") +
xlab("Maternal enviornment")## Warning: Width not defined. Set with `position_dodge(width = ?)`## Warning: Removed 57 rows containing missing values (geom_segment).
#plasticity figure
ggplot(
subset(
LSmeans_tri,
treatment == "Control"),
aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
theme_classic(base_size = 22) +
#geom_point() +
stat_summary(fun=mean, position = "dodge") +
#stat_summary(
# geom="errorbar",
# fun.data= mean_cl_boot,
# width = 0.1, size = 0.2, col = "grey57"
# ) +
# Lines by species using grouping
#scale_color_viridis(discrete = FALSE) +
#scale_fill_manual(values = cbp2)+
stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
ylab("Trichome Density") +
xlab("Maternal enviornment")## Warning: Removed 1 rows containing non-finite values (stat_summary).## Warning: Removed 1 rows containing non-finite values (stat_summary).## Warning: Width not defined. Set with `position_dodge(width = ?)`## Warning: Removed 56 rows containing missing values (geom_segment).
Fitness models: Fitness ~ garden x treatment x elev {.tabset)
Probability of reproduction
Individual
mod_flower<- glmer(flowered~ elev*treatment*mat_treat+ S_init_size + (1|block)+(1|genotype), data = gh2, control=glmerControl(optimizer="optimx", tolPwrss=1e-3, optCtrl=list(method="nlminb")), family=binomial)## Loading required namespace: optimxAnova(mod_flower) # all seperate factors are sig, sig interaction btw mat_treat and elev## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: flowered
## Chisq Df Pr(>Chisq)
## elev 1.6336 1 0.201202
## treatment 10.2749 1 0.001348 **
## mat_treat 10.1867 2 0.006138 **
## S_init_size 71.0068 1 < 2.2e-16 ***
## elev:treatment 1.3670 1 0.242324
## elev:mat_treat 5.4014 2 0.067158 .
## treatment:mat_treat 0.1544 2 0.925696
## elev:treatment:mat_treat 4.8511 2 0.088428 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot (mod_flower) #plot is strange, not sure if ok fit
flower <-predictorEffect("elev", partial.residuals=TRUE, mod_flower)## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictor S_init_size is a one-column matrix that was converted to a vectorplot(flower, lwd=2,xlab="Source Elevation (KM)", ylab="Probability of reproduction", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"),
partial.residuals=list(smooth=TRUE, pch=19, col="black"), ylim=c(0,1))## Warning in plot.eff(flower, lwd = 2, xlab = "Source Elevation (KM)", ylab
## = "Probability of reproduction", : partial residuals are not displayed in a
## multiline plot
#no parents
mod_pr1<- glmer(flowered~ elev*treatment*mat_treat+ S_init_size+(1|block)+(1|genotype), data =subset(gh2, mat_treat!="Parent"), control=glmerControl(optimizer="optimx", tolPwrss=1e-3,optCtrl=list(method="nlminb")), family=binomial)
Anova(mod_pr1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: flowered
## Chisq Df Pr(>Chisq)
## elev 2.7394 1 0.097903 .
## treatment 10.8137 1 0.001008 **
## mat_treat 0.0000 1 0.999811
## S_init_size 47.3706 1 5.876e-12 ***
## elev:treatment 0.0731 1 0.786891
## elev:mat_treat 5.1058 1 0.023846 *
## treatment:mat_treat 0.0705 1 0.790671
## elev:treatment:mat_treat 2.2386 1 0.134606
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1flowerP <-predictorEffect("elev", partial.residuals=TRUE, mod_pr1)## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictor S_init_size is a one-column matrix that was converted to a vectorplot(flowerP, lwd=2,xlab="Source Elevation (KM)", ylab="Probability of reproduction", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"), partial.residuals=list(smooth=FALSE, pch=19, col="black"), ylim=c(0,1))## Warning in plot.eff(flowerP, lwd = 2, xlab = "Source Elevation (KM)", ylab
## = "Probability of reproduction", : partial residuals are not displayed in a
## multiline plot
family level (not working)
Probability of survival
individual
mod_surv<- glmer(survival~elev*treatment*mat_treat + (1|block)+(1|genotype), data = gh2, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e7)), family=binomial)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?Anova(mod_surv) #only treatment is sig## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: survival
## Chisq Df Pr(>Chisq)
## elev 0.3383 1 0.56083
## treatment 9.2003 1 0.00242 **
## mat_treat 1.8540 2 0.39574
## elev:treatment 0.4366 1 0.50878
## elev:mat_treat 0.5099 2 0.77497
## treatment:mat_treat 0.0307 2 0.98477
## elev:treatment:mat_treat 0.0002 2 0.99992
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mod_surv)
survived <-predictorEffect("elev", partial.residuals=FALSE, mod_surv)
plot(survived, lwd=2,xlab="Elevation of origin", ylab="Fitness (survival)", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"),
partial.residuals=list(smooth=FALSE, pch=19, col="black"), ylim=c(0,1))
family level not working
Fecundity
individual
library(glmmTMB)## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.7.21
## Current TMB version is 1.7.22
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)##
## Attaching package: 'glmmTMB'## The following object is masked from 'package:gamlss':
##
## refithurdle <- glmmTMB (Mature_silique_number ~ treatment*mat_treat*elev+S_init_size + (1| block)+(1|genotype), zi=~ treatment*mat_treat*elev+S_init_size + (1| block)+(1|genotype),data = gh2 ,family=truncated_nbinom2)
diagnose(hurdle)## Unusually large coefficients (|x|>10):
##
## zi~(Intercept)
## -12.75103
## zi~mat_treatherb
## 11.52209
## zi~mat_treatParent
## 10.04355
## zi~treatmentHerbivorized:mat_treatherb
## -11.35602
## zi~treatmentHerbivorized:mat_treatParent
## -11.80556
##
## Large negative coefficients in zi (log-odds of zero-inflation),
## dispersion, or random effects (log-standard deviations) suggest
## unnecessary components (converging to zero on the constrained scale);
## large negative and/or positive components in binomial or Poisson
## conditional parameters suggest (quasi-)complete separation
##
##
## Unusually large Z-statistics (|x|>5):
##
## treatmentHerbivorized
## -5.370225
## treatmentHerbivorized:mat_treatParent
## -7.215594
## treatmentHerbivorized:elev
## 24.781566
## treatmentHerbivorized:mat_treatParent:elev
## 11.189258
## zi~treatmentHerbivorized:elev
## -5.921320
## theta_1|block.1
## -6.004913
##
## Large Z-statistics (estimate/std err) suggest a failure of the Wald
## approximation - often also associated with parameters that are at or
## near the edge of their range (e.g. random-effects standard deviations
## approaching 0). While the Wald p-values and standard errors listed in
## summary() are unreliable, profile confidence intervals (see
## ?confint.glmmTMB) and likelihood ratio test p-values derived by
## comparing models (e.g. ?drop1) may still be OK. (Note that the LRT is
## conservative when the null value is on the boundary, e.g. a variance or
## zero-inflation value of 0 (Self and Liang 1987; Stram and Lee 1994;
## Goldman and Whelan 2000); in simple cases the p-value is approximately
## twice as large as it should be.)Anova(hurdle,type="III", component="zi")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mature_silique_number
## Chisq Df Pr(>Chisq)
## (Intercept) 8.4004 1 0.003751 **
## treatment 0.2982 1 0.585032
## mat_treat 8.3993 2 0.015001 *
## elev 9.2922 1 0.002301 **
## S_init_size 56.2618 1 6.344e-14 ***
## treatment:mat_treat 2.9516 2 0.228596
## treatment:elev 0.0285 1 0.865890
## mat_treat:elev 8.4809 2 0.014401 *
## treatment:mat_treat:elev 3.0445 2 0.218217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(hurdle) #zero inflated model is significant, but not sure if i am plotting it?## Family: truncated_nbinom2 ( log )
## Formula:
## Mature_silique_number ~ treatment * mat_treat * elev + S_init_size +
## (1 | block) + (1 | genotype)
## Zero inflation:
## ~treatment * mat_treat * elev + S_init_size + (1 | block) + (1 | genotype)
## Data: gh2
##
## AIC BIC logLik deviance df.resid
## 1961.9 2108.2 -949.9 1899.9 799
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## block (Intercept) 0.03159 0.1777
## genotype (Intercept) 0.21906 0.4680
## Number of obs: 830, groups: block, 12; genotype, 19
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## block (Intercept) 0.9225 0.9605
## genotype (Intercept) 1.1907 1.0912
## Number of obs: 830, groups: block, 12; genotype, 19
##
## Dispersion parameter for truncated_nbinom2 family (): 4.47
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.50261 1.77175 3.106 0.0019
## treatmentHerbivorized -0.47436 2.54744 -0.186 0.8523
## mat_treatherb -1.83915 1.48197 -1.241 0.2146
## mat_treatParent -0.84673 1.60535 -0.527 0.5979
## elev -1.09470 0.58242 -1.880 0.0602
## S_init_size 0.09895 0.07036 1.406 0.1597
## treatmentHerbivorized:mat_treatherb 1.45182 3.64446 0.398 0.6904
## treatmentHerbivorized:mat_treatParent -0.50276 3.62770 -0.139 0.8898
## treatmentHerbivorized:elev 0.03551 0.88005 0.040 0.9678
## mat_treatherb:elev 0.64279 0.49719 1.293 0.1961
## mat_treatParent:elev 0.31005 0.53393 0.581 0.5615
## treatmentHerbivorized:mat_treatherb:elev -0.47540 1.26443 -0.376 0.7069
## treatmentHerbivorized:mat_treatParent:elev 0.11147 1.24729 0.089 0.9288
##
## (Intercept) **
## treatmentHerbivorized
## mat_treatherb
## mat_treatParent
## elev .
## S_init_size
## treatmentHerbivorized:mat_treatherb
## treatmentHerbivorized:mat_treatParent
## treatmentHerbivorized:elev
## mat_treatherb:elev
## mat_treatParent:elev
## treatmentHerbivorized:mat_treatherb:elev
## treatmentHerbivorized:mat_treatParent:elev
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.7510 4.3994 -2.898 0.00375
## treatmentHerbivorized 3.0553 5.5953 0.546 0.58503
## mat_treatherb 11.5221 4.1898 2.750 0.00596
## mat_treatParent 10.0436 4.2318 2.373 0.01763
## elev 4.3293 1.4202 3.048 0.00230
## S_init_size -1.2138 0.1618 -7.501 6.34e-14
## treatmentHerbivorized:mat_treatherb -11.3560 7.6944 -1.476 0.13997
## treatmentHerbivorized:mat_treatParent -11.8056 7.8773 -1.499 0.13395
## treatmentHerbivorized:elev -0.3140 1.8594 -0.169 0.86589
## mat_treatherb:elev -3.8702 1.3659 -2.833 0.00461
## mat_treatParent:elev -3.0781 1.3756 -2.238 0.02524
## treatmentHerbivorized:mat_treatherb:elev 3.9197 2.5866 1.515 0.12967
## treatmentHerbivorized:mat_treatParent:elev 3.9602 2.6349 1.503 0.13284
##
## (Intercept) **
## treatmentHerbivorized
## mat_treatherb **
## mat_treatParent *
## elev **
## S_init_size ***
## treatmentHerbivorized:mat_treatherb
## treatmentHerbivorized:mat_treatParent
## treatmentHerbivorized:elev
## mat_treatherb:elev **
## mat_treatParent:elev *
## treatmentHerbivorized:mat_treatherb:elev
## treatmentHerbivorized:mat_treatParent:elev
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fruit <-predictorEffect("elev", partial.residuals=TRUE, hurdle,component="zi")## Warning in Effect.glmmTMB(ans, mod, x.var = 1, xlevels = xlevels, ...):
## overriding variance function for effects: computed variances may be incorrect## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictor S_init_size is a one-column matrix that was converted to a vectorplot(fruit, lwd=2,xlab="Source Elevation (KM)", ylab="Fruit", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"), partial.residuals=list(smooth=FALSE, pch=19, col="black"), ylim=c(0,25))## Warning in plot.eff(fruit, lwd = 2, xlab = "Source Elevation (KM)", ylab =
## "Fruit", : partial residuals are not displayed in a multiline plot
Family level
LAR
individual
#have to curate the dataset to focus on herb treatment and remove NAs
#LAR 3 has the most damage, LAR 13 is the most recent
LAR_data <- subset(gh2, treatment == "Herbivorized")
LAR_3 <- drop_na(LAR_data,LAR_prop3) #this removes the rows without LAR values
LAR_3 <- dplyr::select(LAR_3, LAR_prop3, elev, genotype, treatment, mat_treat, exp_ID, S_init_size, block)
gamlss_mod1<- gamlss (formula=LAR_prop3~elev*mat_treat+ S_init_size+ random(block)+random(genotype), family=BE, data= LAR_3)## GAMLSS-RS iteration 1: Global Deviance = -317.6556
## GAMLSS-RS iteration 2: Global Deviance = -333.6844
## GAMLSS-RS iteration 3: Global Deviance = -334.1667
## GAMLSS-RS iteration 4: Global Deviance = -334.1874
## GAMLSS-RS iteration 5: Global Deviance = -334.1902
## GAMLSS-RS iteration 6: Global Deviance = -334.191summary(gamlss_mod1)## ******************************************************************
## Family: c("BE", "Beta")
##
## Call: gamlss(formula = LAR_prop3 ~ elev * mat_treat + S_init_size +
## random(block) + random(genotype), family = BE, data = LAR_3)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: logit
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.270656 0.774938 -1.640 0.10192
## elev 0.326973 0.249718 1.309 0.19122
## mat_treatherb -0.038300 1.110804 -0.034 0.97251
## mat_treatParent 0.066532 1.096706 0.061 0.95166
## S_init_size -0.124692 0.039976 -3.119 0.00196 **
## elev:mat_treatherb -0.024684 0.358378 -0.069 0.94512
## elev:mat_treatParent -0.004827 0.353084 -0.014 0.98910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: logit
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.58078 0.04473 -12.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 399
## Degrees of Freedom for the fit: 29.87062
## Residual Deg. of Freedom: 369.1294
## at cycle: 6
##
## Global Deviance: -334.191
## AIC: -274.4498
## SBC: -155.297
## ******************************************************************plot(gamlss_mod1)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.006089556
## variance = 0.9862509
## coef. of skewness = 0.2791661
## coef. of kurtosis = 5.366551
## Filliben correlation coefficient = 0.9871931
## ******************************************************************LAR_13 <- drop_na(LAR_data,LAR_prop13) #this removes the rows without LAR values
LAR_13 <- dplyr::select(LAR_13, LAR_prop13, elev, genotype, treatment, mat_treat, exp_ID, S_init_size, block)
gamlss_mod2<- gamlss (formula=LAR_prop13~elev*mat_treat+ random(block)+random(genotype), family=BEZI, data= LAR_13)## GAMLSS-RS iteration 1: Global Deviance = -608.1084
## GAMLSS-RS iteration 2: Global Deviance = -745.3908
## GAMLSS-RS iteration 3: Global Deviance = -864.3819
## GAMLSS-RS iteration 4: Global Deviance = -963.2275
## GAMLSS-RS iteration 5: Global Deviance = -1025.301
## GAMLSS-RS iteration 6: Global Deviance = -1049.247
## GAMLSS-RS iteration 7: Global Deviance = -1054.845
## GAMLSS-RS iteration 8: Global Deviance = -1055.889
## GAMLSS-RS iteration 9: Global Deviance = -1056.105
## GAMLSS-RS iteration 10: Global Deviance = -1056.16
## GAMLSS-RS iteration 11: Global Deviance = -1056.175
## GAMLSS-RS iteration 12: Global Deviance = -1056.18
## GAMLSS-RS iteration 13: Global Deviance = -1056.181
## GAMLSS-RS iteration 14: Global Deviance = -1056.182summary(gamlss_mod2)## ******************************************************************
## Family: c("BEZI", "Zero Inflated Beta")
##
## Call: gamlss(formula = LAR_prop13 ~ elev * mat_treat + random(block) +
## random(genotype), family = BEZI, data = LAR_13)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: logit
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.6303 0.9162 -2.871 0.0044 **
## elev -0.1151 0.2955 -0.390 0.6972
## mat_treatherb -0.9045 1.2696 -0.712 0.4768
## mat_treatParent 0.7480 1.2484 0.599 0.5496
## elev:mat_treatherb 0.2692 0.4091 0.658 0.5111
## elev:mat_treatParent -0.2596 0.4023 -0.645 0.5193
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.3235 0.0894 37.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: logit
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.6009 0.2262 -11.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 304
## Degrees of Freedom for the fit: 22.51564
## Residual Deg. of Freedom: 281.4844
## at cycle: 14
##
## Global Deviance: -1056.182
## AIC: -1011.151
## SBC: -927.4594
## ******************************************************************plot(gamlss_mod2)
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.196609
## variance = 0.6954236
## coef. of skewness = -0.1419971
## coef. of kurtosis = 2.991946
## Filliben correlation coefficient = 0.9931949
## ******************************************************************visreg(gamlss_mod1, 'elev', by= "mat_treat", overlay = FALSE, type="conditional",
#scale = "response",
xlab="Elevation (m)", ylab="Percentage leaf area herbivorized", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col="blue"),
line=list(col=grey(c(0,0.8))),
points=list(cex=1.5,col=grey(c(0,0.8)))) ## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
visreg(gamlss_mod2, 'elev', by= "mat_treat", overlay = FALSE, type="conditional",
#scale = "response",
xlab="Elevation (m)", ylab="Percentage leaf area herbivorized", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
fill=list(col="blue"),
line=list(col=grey(c(0,0.8))),
points=list(cex=1.5,col=grey(c(0,0.8)))) ## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## New way of prediction in random() (starting from GAMLSS version 5.0-6)